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Computation of Pressurized Gas Bearings Using CE/SE Method 

Sorin Cioc, Florin Dimofte, and Theo G. Keith, Jr. 

University of Toledo 
Toledo, Ohio 43606 

David P. Fleming 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

The space-time conservation element and solution element (CE/SE) method is extended to 
compute compressible viscous flows in pressurized thin fluid films. This numerical scheme has 
previously been used successfully to solve a wide variety of compressible flow problems, including 
flows with large and small discontinuities. In this paper, the method is applied to calculate the pressure 
distribution in a hybrid gas journal bearing. The formulation of the problem is presented, including the 
modeling of the feeding system. The numerical results obtained are compared with experimental data. 
Good agreement between the computed results and the test data were obtained, and thus validate the 
CE/SE method to solve such problems. 

NOMENCLATURE 

a Radius of feeding orifice 

a F ,b F ,c F Coefficients in Taylor series expression of function / 

a Cr ,b G ,c G Coefficients in Taylor series expression of function g 

C Radial clearance 

D Diameter of feeding pocket 

e Eccentricity 

F fi + gk Vector flux term 

/ ircumferential flux term, see Eq. (5) 

g Axial flux term, see Eq. (5) 

h Film thickness (dimensional) 

h h/C Non-dimensional film thickness 

i ,k unit vectors in circumferential and in axial directions, respectively 

L Length of bearing 

n Adiabatic coefficient 

n w Number of waves (for wave bearings) 

ns Number of feeding holes 

p Fluid pressure 

p pj p 0 Non-dimensional pressure 

p 0 Atmospheric (reference) pressure 

p s Supply pressure 

0 Mass flow rate through one feeding hole 

O Non-dimensional mass flow rate through one feeding hole (see Eq. (17)) 

R Journal bearing radius 
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h p Dependent variable in governing equation 
co R Velocity in circumferential direction 
pco f R 
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C 


Non-dimensional circumferential velocity 


Circumferential coordinate 

x/(2nR) Non-dimensional circumferential coordinate 

Position of maximum fdm thickness, measured in negative direction of axis x 


Position of wave peak (for wave bearings), measured in negative direction of axis x 
Axial coordinate 

z/(2nR) Non-dimensional axial coordinate 
Isentropic exponent; y = 1.405 for air 
Grid spacing in longitudinal direction 
Time step 

Weight parameter characterizing one form of artificial dissipation, or 

e/C Eccentricity ratio 

Wave amplitude (for wave bearings) 

Fluid viscosity 
Fluid density 

p/p 0 Non-dimensional fluid density 

Atmospheric (reference) density 
Angular velocity of the journal bearing 


Subscripts and superscripts, other than shown above 

( ) 0 Value of the variable at the point O 

( ) ; Supply pocket index 

( )x ’( ) z >( ) f Partial derivative with respect to x, z, and t, respectively 
( )” Time step 


INTRODUCTION 

Gas lubrication, though preceded by occasional experimental work since the mid 19 th century, 
experienced a strong development in the “Golden Era” of gas lubrication, which started in the last years 
of World War II, and ended in the first part of the 1970’s (Pan, 1990). In this period much important 
work was published: Fuller (1956), Lund (1964, 1967), Castelli and Pirvics (1968), Constantinescu 
(1969). More recently, significant advancements were reported in gas film modeling including 
rarefaction effects (Wu and Bogy, 2001), in numerical methods applicable to high-speed bearings (Faria 
and San Andres, 2000), and particularly in the treatment of complex geometries, including 
discontinuities (Bonneau, Huitric and Toumerie, 1993). 
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The space-time conservation element and solution element (CE/SE) method was proposed by 
Chang and To (1991). Over the past several years it has been utilized in a number of fluid flow 
applications that involve shock waves, contact discontinuities, acoustic waves, vortices and chemical 
reactions. One of its main features is that it can simultaneously capture small and large discontinuities 
(such as sound waves and shock waves) without introducing numerical oscillations into the solution, as 
shown by Chang, Wang, and Chow (1999), and Chang et al. (1998). 

Compressible viscous flow in pressurized thin fluid films, with application in hybrid gas bearings, 
can encounter large pressure gradients due to the feeding system or to the large peripheral velocities of 
the bearing, hi these conditions, the computational methods based on standard finite difference methods 
or classic finite volume methods may prove inadequate, having convergence problems or inducing 
oscillations into the solution. Thus, a method that is conceptually simple, is second order accurate for 
the entire domain, and is able to naturally deal with large gradients and/or discontinuities in the solutions 
without introducing numerical oscillations or smearing, is welcome. 

ANALYSIS 

The two-dimensional, transient, Reynolds equation, written for a Newtonian compressible fluid in 
laminar flow is, 


dp h d f phU 
dt 2 


ph' dp' + d ' 

12 pdxj 


ph 3 d p N 

12 ( 0 . dz J 


( 1 ) 


The flow is considered polytropic, i.e., 

~ = const . , (2) 

P 

where the polytropic exponent n = 1 for isothermal flow, n = 1.405 for adiabatic flow, or can have other 
values for general polytropic flows. 

A more suitable form of the Reynolds equation for numerical formulation is obtained using a 
new variable u that is the product of the non-dimensional film thickness and the pressure, i.e., 

u = hp (3) 


In terms of u, in non-dimensional variables, the Reynolds equation can be written as 


du df dg 
— + -^~ + — 2 - = 0 , 
dt dx dz 


( 4 ) 


where the flux terms / and g are, 
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( 5 ) 


All partial derivatives considered in Eqs. (4) and (5) are carried out relative to non-dimensional variables 
t,x,z. In the following, in order to simplify the expressions, the non-dimensional notation (upper bar) 
will be dropped, and all variables will be implicitly considered in non-dimensional form. 

Consider a triangular mesh that covers the (x,z) spatial domain. One triangle BCD and its three 
neighboring elements are shown in Fig. 1. Point A is the centroid of the triangle BCD, while points E, F 
and G are the centroids of the neighboring triangles BCH, CDI and Bill, respectively. The CE/SE 

i 

n- 1 — 

method calculates the values of the dependent variables u,u x ,u : for point A at the time step t = t 2 

using the corresponding values of the same variables for the points E, F and G at the time step t = t n . In 
order to calculate the three unknowns at the new time step, a system of three equations will be derived. 

Consider the quadrilateral ABEC. Simultaneously integrating Eq. (4) over the surface of this 

i 

n+— 

quadrilateral and in time, between time steps t n and t 2 (see Fig. 2), yields 


ABEC t n 


II 1^-1 1/ 


f ABEC 


df , dg 

dx dz 


dad? = 0 . 


( 6 ) 


Performing the time integration for the first term and transforming the surface integration into a contour 
integration for the second term (using Green’s theorem) produces 


If 

ABEC 



i 

da+ j j>F-/7dsd/ = 0, 


t" ABEC 


(7) 


where n is the unit vector normal to the contour, oriented outwards and 

F = fi+gk (8) 

is the vector in the (x, z) plane characterized by the Cartesian unit vectors (/ .k). Functions / and g are 
given by Eq. (5). Equation (7) implies conservation of flux in the three-dimensional space (x,z,t). 
Functions u,f g are next written with linear approximations using first order Taylor expansions, i.e., 

u = Uq + (u x ) 0 (x - x 0 ) + {it z \ (z-z 0 )+ (u t (9) 


/ = /o + 


f df} 

\ du lo 




f df? 

K du Do 


k -(«,)<>] 


+ 


f dj^ 

\ du D 0 


k-k)o ]= a F u+b F u x +c F , 


( 10 ) 
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Substituting Eq. (9) into Eqs. (10) and (11) yields linear expressions for / and g as functions of 
(x,z,t). Li Eqs. (10) and (11), coefficients a F ,b F ,c F ,a G ,b G ,c G are considered constant when 
integrating Eq. (7), and are known as functions of u 0 , (u x ) 0 , (u z ) 0 . In Eq. (9), the time derivative can be 
calculated as function of the space derivatives using Eq. (3), i.e., 
k )o = “(A- )o - igz )o = ~ a F k )o - a G O', )o • ( ■ 1 2 ) 

f ^ 

Equations (8-12) are then substituted into Eq. (7). Point x A .,z A ,,t 2 is used as the Taylor 


expansion point for the expressions of u 2 ,/ and g on the contour segments AB and CA, while point 
(x E ,,z E ,,t n ) is used as the Taylor expansion point for the expressions of //",/ and g on the contour 
segments BE and EC. Thus, a first equation with three unknowns, the values u,u x ,u z at the new half 

f 

time step x A .,z A .,t 2 , is obtained. The coordinates of points A' and E' are selected in a suitable 


way, as will be shown later. The equation is linear and has the general form 

i i i 

n+— y v.«H y N.W+— y \ y \ 

a x u A , 2 +b x (u x ) A , 2 + c x {u : ) A ' 2 + d i u E'+ e MX + fi( u X + gi = °- ( 13a ) 

Two other similar equations are obtained using the same procedure for the conservation elements 
ACFD and ADGB. These equations have the general form 


i i i 

n+— y \ n+ ~ / \ n+ ~ ( \n ( \n 

a 2 U A' 2 + b 2 (».v )jt 2 + C 2 (», h 2 + d 2 u” F . + e 2 (u x ) F . + f 2 {u 2 ) F . + g 2 = 0 


(13b) 


1 1 1 

n+— y \ n+ ~ / \ n+ ~ ( \n ( \n 

a,u A , 2 + b } (u x ) A , 2 + c 3 (u 2 ),, 2 + d 3 u n G , + e 3 (u x ) G , + f 3 (u 2 ) a ,+g 3 =0. (13c) 

The linearized system formed by Eqs. (13 a), (13b), and (13 c) can be solved using an iterative 

i i i 

method; note that the coefficients a i ,b i ,c i ,i = 1,2,3 are functions of the unknowns u A , 2 , \u x ) A , 2 ,{u 2 ) A , 2 . 
An alternate approach is to choose the Taylor series expansion points' as the center of the hexagon 
BEC.FDG. The other three Taylor series expansion points E',F',G' can be chosen arbitrarily; however, 
in order to maintain consistency, they are selected as the centers of the corresponding hexagons formed 
around the neighboring triangular elements. 

Adding Eqs. (13a), (13b), and (13c) yields a new equation that represents the flux conservation 
over the hexagon and over one half time step. When point A' is the centroid of the hexagon BECFDG, 
this equation has a simpler form given by 
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(14) 


1 

^sum^A' ^ ~^-d]Ugr + £q (“, )f + f ( U z Ye + ^2 W F' + *2 (“, )f' + fl ( U z Yf' 

d Aa + e 3 (m x )”, + f 3 (u z )”, + g,„ m = 0 , 


where a ram is given by 


Cl sum A abeC + -A ACFD + ^ADGB _ ^ BECFDG • 


(15) 


i 

n+— 

Equation (14) has only one unknown, u A , 2 , and can easily be solved explicitly. It is also 
important to note that all the coefficients in Eq. (14) depend only on the geometry (coordinates of the 
points) and the values of the dependent variables at the previous half time step so that an iterative 

i 

— 

method is not needed. After calculating the value u A , 2 , the values of the other two dependent variables 

1 1 

/ \«H — / \W+— 

[u x ) A , 2 and \u z ) A , 2 can be calculated using any two of the Eqs. (13a), (13b), and (13c). This is also 
called the a scheme. 


Note that the a scheme introduces “no significant damping” (Chang, Wang, & Chow, 1999; 
Chang et al., 1998), so that some form of artificial dissipation is generally necessary. The scheme can be 
simplified and simultaneously stabilized by calculating the space derivatives in a different way. The 
scheme thus obtained is called the a- s - a - p scheme. In this scheme, the derivatives are calculated as 
weighted averages between the derivatives calculated from the governing equations, as shown above 
(the a scheme), the derivatives calculated using 2-D central difference finite difference formulae (weight 
parameter s) and the derivatives using 2-D side finite differencing (weight parameter P). Parameter a is 
the power index used in the computation of the non-linear weighted average that employs 2-D side finite 
differencing. The complete formulation of the way the derivatives are calculated in the a-s-a-p 
scheme can be found in Chang et al. (1998). It is important that, for a certain value of the weighting 
parameter s (s = 0.5), the a-s-a-p scheme eliminates the necessity of calculating the space 
derivatives from the governing equations, thus the method becomes purely explicit. 

The feed system, formed by a number ns of orifices with the general geometry shown in Fig. 3, 
can be modeled using a generally accepted formula (Lund, 1964, 1967), that links the mass flow ratio to 
the feed system geometry and the pressures at the ends of the feed system considering both effects of the 
orifice restrictor and the inherent restrictor, i.e.. 


o 




f 2 y 

' a ' 


dh 


o. 


(16) 


where the non-dimensional mass flow 0 is a fraction C D of the ideal mass flow, 
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When no restricting orifice is present (i.e., an inherent restrictor), Eq. (16) becomes 

Q = TidhyJ p s p s Q (18) 

The feeding system is introduced into the bearing computation through the boundary conditions. 
Thus, on each feeding pocket contour, the mass flow rates calculated from the bearing equations and 
from the feeding system, respectively, must be equal 

fe \earing ~ (— i ) feeding ’ ? — b 2, •••,/« . (19) 

Equation (19) represents a nonlinear set of ns equations, where the pressures in the bearing at the 
feeding holes are the unknowns p, ,i = \,2,---,ns . This system is solved iteratively using Newton’s 
method. 

RESULTS AND DISCUSSION 

A computer code has been developed based on the described method, using the a-s-a-p 
scheme with s = 0.5 . This code has been tested for some simple cases where experimental data were 
available, considering isothermal flow (n = l) . 

The first case considered was the circular gas bearing without any feeding system. The results are 
shown in Fig. 4, where the same bearing with aspect ratio L/(2R) = 2.0 is examined for two bearing 
numbers A = (6p©A")/(p 0 C 2 ) and for two eccentricities. There are no significant differences between 
the experimental and calculated pressure distributions, at least in the middle plane where the 
experimental results are available. The relative differences between the experimental and calculated 
non-dimensional loads C, = Load /(p 0 2LR) are 3.5% and 4.0%; both represent improvements over the 
theoretical results shown by Constantinescu (1969). 

The second case considered is the wave journal bearing (Dimofte, 1995) without any feed system. 
The fluid film thickness in an aligned wave bearing is given as 

h =l + scos27i(x + x 0 )+8 w cos[27i:» M , (x + xj. (20) 

In the last term on the right hand side of Eq. (19), e w is the wave amplitude, n w is the number of waves 
and x w is the wave origin relative to the origin of the x axis. For the present case, the bearing diameter 
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and bearing length are 50 mm, the radial clearance is 0.02 mm, n w = 3, s w = 0.3, and the relative 
minimum film thickness is h mm = 0.3 for all cases. The bearing number A is 3.566. 

The results are compared with a finite difference based code built by Dimofte (1995). Different 
bearing positions ( x w ) have been tested and the results (calculated load magnitude and position) are 
presented in Table 1. The results show very good agreement between the two methods. 

The third case considered is a pressurized gas bearing without rotation and with zero eccentricity. 
Details of the bearing geometry and working conditions are presented in Table 2. Figures 5a and 5b 
show the calculated and experimental pressure distributions for two values of the radial clearance, 
corresponding to the subsonic flow regime (obtained when C = 1 2.7 \im ), and the choked flow regime 
(obtained whenC = 3 1.75 pm ) in the feed system, in two longitudinal planes situated at the jet position 
(Fig. 5a) and at half distance between jets (Fig. 5b). Pressure peaks are visible at the jet positions; also 
nearly constant pressure is obtained between the supply planes. The predicted pressure peaks at the 
position of jet are higher than the experimental values (this difference is more visible for subsonic inlet 
flow, however it is present in both cases). This is due to the difficulty of measuring the local pressure at 
the feeding orifice position. The value for the correction factor Co is 0.8 for subsonic inlet flow, and 
0.85 for choked flow. 

Finally, the code was applied to different configurations where experimental or computed data 
were not available. 

Figure 6 shows the pressure distribution (sections at the jet and at the mid-jet positions) obtained 
for five supply planes for the same bearing as in the previous case for C = 12.7 pm . The supply planes 
are equally distanced with respect to each other and with the bearing ends. It is seen that, although the 
external pressure is the same for all supply holes, the pressure that develops in the pockets (inside the 
bearing) is not the same for all supply planes. The pressure distribution between two consecutive feeding 
planes has an almost linear form. Furthermore, at the central supply plane, the peak pressures are not as 
prominent as the peak pressures at the other supply planes. This suggests that the central supply plane 
does not have an important contribution for the general pressure distribution inside the bearing. 

Figure 7 shows the pressure distribution in a half bearing for the same bearing geometry as in the 
previous two cases, but in a running condition (40,000 rpm, A = 25.76), with a relative eccentricity 
s = 0.5 . The value C. D =0.8 was used. The pressure distribution is very complex. It can be seen that for 
many feeding holes the flow is inverted, i.e., is directed from the bearing towards the feeding system, 
because the pressure inside the bearing is higher than the supply pressure. 

CONCLUSIONS 

The CE/SE computational method has been extended for the first time to calculate compressible 
viscous flow in pressurized thin fluid films with restrictors in series. Formulation of the method was 
presented along with numerical results obtained for both non-pressunzed and pressurized bearings. The 
results were compared with existing expenmental and computed data. The results demonstrate the 
ability of the method to accurately predict pressure distribution in such flows. 

While most numerical methods handle space and time differencing terms in the governing 
equations separately, this relatively new scheme treats them in a unified way. This results in very good 
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accuracy even though the unknowns are considered locally linear and even when a relatively course grid 
is used. Also, the entire structure of the method, as well as the developed code, is relatively simple. 
Unlike most numerical schemes, no special treatment is necessary for discontinuities, providing that the 
governing equations are written in strong conservative form. However, the method is limited to time 
dependent equations. Steady state results, such as those discussed in this work, can be obtained only 
after the stabilization of the unsteady solution starting from initial conditions. 

Future work will focus on predicting dynamic characteristics of pressurized gas bearings. 
Improving the feeding system with restrictors in series modeling will also be considered, as well as 
improving the computation time. Cases including discontinuities will also be considered. 
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Table 1. - Comparison between calculated data with CE/SE method and with finite 

difference method 


Eccentricity 

«w 

(deg) 

Load CE/SE 
(N) 

Load FD 
(N) 

Load angle 
CE/SE (deg) 

Load angle 
FD (deg) 

0.502 

40 

177 

177 

40.49 

40.18 

0.460 

32 

170 


39.09 

38.8 

0.433 

24 

166 

165 

36.37 

36.11 

0.416 

16 

163 

162 

32.94 

32.69 

0.404 

8 

158 

158 

29.30 

29.06 

0.400 

0 

154 

154 

25.61 

25.38 

0.404 

-8 

151 

151 

22.08 

21.86 

0.416 

-16 

149 

148 

18.95 

18.75 


Table 2. - Pressurized gas bearing geometry and working conditions 


Bearing length (mm) 

117.5 

Bearing diameter (mm) 

60.4 

Supply planes 

2 

Supply plane position (mm) 

12.7 

Holes/supply plane 

14 

Orifice diameter (mm) 

0.16 

Pocket diameter (mm) 

0.9 

Supply pressure (Pa) 

5.514xl0 3 

Injection angle (deg) 

90 
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Fig. 1 - Triangular mesh element and its neighbors. 



x 

Fig. 2 - Conservation volume in ( x,z,t ) space. 
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A = 0.45,s exp =0.41,Cexp =O.1496,0 exp =65° 

A = 1.569, 8 exp =0.51, C exp =0.545, 0 exp =37.5° 

«W=0.40,^= 0.1548,0^=69° 

8^ =0.51,^ =0.567,0^ =38.0° 


Fig. 4 - Comparison between numerical and experimental results for pressure distribution 

in gas journal bearing. 
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Fig. 6 - Pressure distribution for circular bearing without eccentricity; five supply planes. 



Fig. 7 - Pressure distribution for circular bearing with eccentricity; five supply planes. 
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